function HPs = dropout_update_HP_forward(H, HPs, Ks, beta, gamma, eta,  mid_layer_clusters, mid_layer_width, layer_connections, layer_connections_backward)
nLayer = length(mid_layer_clusters);
[nSmp, ~, ~] = size(Ks);
%**********************************************************************
% update H1_1->H1_2->H1_m->H2_1->H2_2...->H2_m->...->H*
%**********************************************************************

if nLayer == 1
    iLayer = 1;
    layer_pre_forward = layer_connections{iLayer};
    gamma_pre = gamma{1, iLayer};
    eta_a = eta{1, iLayer};
    
    tmpHPt = cell(mid_layer_width(iLayer), 1);
    
    for i1 = 1:mid_layer_width(iLayer)
        tmp1 = zeros(nSmp);
        current_input_output = layer_pre_forward{i1, 1};
        gamma_input_output = gamma_pre{i1, 1};
        for i2 = 1:length(current_input_output)
            tmp1 = tmp1 + gamma_input_output(i2) * Ks(:, :, current_input_output(i2));
        end
        
        tmp2 = beta(i1) * (H * H');
        
        tmp = eta_a(i1) * tmp1 + tmp2 + eps;
        tmp= (tmp + tmp')/2;
        opt.disp = 0;
        [H_a_b,  ~] = eigs(tmp, mid_layer_clusters(iLayer), 'LA', opt);
        tmpHPt{i1, 1} = H_a_b;
    end
    HPs{iLayer} = tmpHPt;
else
    
    for iLayer = 1:nLayer
        if iLayer == 1
            layer_pre_forward = layer_connections{iLayer};
            gamma_pre = gamma{1, iLayer};
            eta_a = eta{1, iLayer};
            
            layer_next_backward = layer_connections_backward{iLayer + 1};
            gamma_next = gamma{1, iLayer + 1};
            eta_b = eta{1, iLayer + 1};
            
            tmpHPt = cell(mid_layer_width(iLayer), 1);
            
            for i1 = 1:mid_layer_width(iLayer)
                tmp1 = zeros(nSmp);
                current_input_output = layer_pre_forward{i1, 1};
                gamma_input_output = gamma_pre{i1, 1};
                for i2 = 1:length(current_input_output)
                    tmp1 = tmp1 + (gamma_input_output(i2) * Ks(:, :, current_input_output(i2)));
                end
                               
                
                tmp2 = zeros(nSmp);
                next_input_output = layer_next_backward{i1, 1};
                HPs_next = HPs{1, iLayer + 1};
                for i2 =  1:length(next_input_output)
                    gamma_input_output = gamma_next{next_input_output(i2)};
                    for i3 = 1 : length(gamma_input_output)
                        if gamma_input_output(i3) ~= -1
                            gamma_input_output_idx = i3;
                            break;
                        else
                            continue;
                        end
                    end
                    gamma_input_output_value = gamma_input_output(gamma_input_output_idx);
                    tmp2 = tmp2 + eta_b(next_input_output(i2)) * gamma_input_output_value * (HPs_next{next_input_output(i2)}  * HPs_next{next_input_output(i2)}');
                    gamma_next{next_input_output(i2)}(gamma_input_output_idx) = -1;
                end
                
               
                tmp = eta_a(i1) * tmp1 + tmp2 + eps;
                tmp= (tmp + tmp')/2;
                opt.disp = 0;
                [H_a_b,  ~] = eigs(tmp, mid_layer_clusters(iLayer), 'LA', opt);
                tmpHPt{i1, 1} = H_a_b;
            end
            
            flag = 0;
            for i1 = 1 : size(gamma_next, 1)
                temp = gamma_next{i1};
                for i2 = 1 : length(temp)
                    if temp(i2) ~= -1
                        flag = 1;
                    end
                end
            end
            assert(flag == 0)
            
        elseif iLayer == nLayer
            
            layer_pre_forward = layer_connections{iLayer};
            gamma_pre = gamma{1, iLayer};
            eta_a = eta{1, iLayer};
            
            
            tmpHPt = cell(mid_layer_width(iLayer), 1);
            
            for i1 = 1:mid_layer_width(iLayer)
                tmp1 = zeros(nSmp);
                current_input_output = layer_pre_forward{i1, 1};
                gamma_input_output = gamma_pre{i1, 1};
                HPs_current = HPs{1, iLayer - 1};
                for i2 = 1:length(current_input_output)
                    tmp1 = tmp1 + gamma_input_output(i2) * (HPs_current{current_input_output(i2), 1} * HPs_current{current_input_output(i2), 1}' );
                end
                
                tmp2 = beta(i1) * (H * H');
                tmp = eta_a(i1) * tmp1 + tmp2 + eps;
                tmp= (tmp + tmp')/2;
                opt.disp = 0;
                [H_a_b,  ~] = eigs(tmp, mid_layer_clusters(iLayer), 'LA', opt);
                tmpHPt{i1, 1} = H_a_b;
            end
        else
            layer_pre_forward = layer_connections{iLayer};
            gamma_pre = gamma{1, iLayer};
            eta_a = eta{1, iLayer};
            
            layer_next_backward = layer_connections_backward{iLayer + 1};
            gamma_next = gamma{1, iLayer + 1};
            eta_b = eta{1, iLayer+1};
            
            tmpHPt = cell(mid_layer_width(iLayer), 1);
            
            for i1 = 1:mid_layer_width(iLayer)
                tmp1 = zeros(nSmp);
                current_input_output = layer_pre_forward{i1, 1};
                gamma_input_output = gamma_pre{i1, 1};
                HPs_current = HPs{1, iLayer + 1};
                for i2 = 1:length(current_input_output)
                    tmp1 = tmp1 + gamma_input_output(i2) * (HPs_current{current_input_output(i2), 1} * HPs_current{current_input_output(i2), 1}' );
                end
                
                tmp2 = zeros(nSmp);
                next_input_output = layer_next_backward{i1, 1};
                HPs_next = HPs{1, iLayer + 1};
                for i2 =  1:length(next_input_output)
                    gamma_input_output = gamma_next{next_input_output(i2)};
                    for i3 = 1 : length(gamma_input_output)
                        if gamma_input_output(i3) ~= -1
                            gamma_input_output_idx = i3;
                            break;
                        else
                            continue;
                        end
                    end
                    gamma_input_output_value = gamma_input_output(gamma_input_output_idx);
                    tmp2 = tmp2 + eta_b(next_input_output(i2)) * gamma_input_output_value * (HPs_next{next_input_output(i2)}  * HPs_next{next_input_output(i2)}');
                    gamma_next{next_input_output(i2)}(gamma_input_output_idx) = -1;
                end
                
                tmp = eta_a(i1) * tmp1 + tmp2 + eps;
                tmp= (tmp + tmp')/2;
                opt.disp = 0;
                [H_a_b,  ~] = eigs(tmp, mid_layer_clusters(iLayer), 'LA', opt);
                tmpHPt{i1, 1} = H_a_b;
            end
        end
        HPs{iLayer} = tmpHPt;
    end
    
end